
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd


def get_dft(x):
    N = len(x)
    t = np.linspace(0, 1, N+1)[0:N]
    n_freqs = int(np.ceil(N/2))
    cos_sums = np.zeros(n_freqs)
    sin_sums = np.zeros(n_freqs)
    for i in range(n_freqs):
        c = np.cos(2*np.pi*i*t)
        s = np.sin(2*np.pi*i*t)
        cos_sums[i] = np.sum(c*x)/(N/2)
        sin_sums[i] = np.sum(s*x)/(N/2)
    return cos_sums, sin_sums


N = 100
n = np.arange(N)
x = 2*np.cos(2*np.pi*(3)*n/N)
x += np.sin(2*np.pi*(5)*(n-2)/N)
c, s = get_dft(x)

plt.figure(figsize=(12, 3))
plt.plot(x)
plt.title("Original Signal")

plt.figure(figsize=(8, 4))
plt.subplot(121)
plt.stem(c[0:10])
plt.title("Cosine Components")
plt.subplot(122)
plt.stem(s[0:10])
plt.title("Sine Components")

Text(0.5, 1.0, 'Sine Components')


def get_idft(c, s):
    """
    Given the cosine and sine components, reconstruct the signal x
    in the time domain
    """
    N = 2*len(c)
    n = np.arange(N)
    x = np.zeros(N)
    n_freqs = len(c)
    for k in range(n_freqs):
        x += c[k]*np.cos(2*np.pi*(k/N)*n)
        x += s[k]*np.sin(2*np.pi*(k/N)*n)
    return x


x_recon = get_idft(c, s)
plt.figure(figsize=(12, 4))
plt.plot(x_recon)
plt.plot(x, linestyle='--')
plt.legend(["Reconstructed Signal", "Original Signal"])

<matplotlib.legend.Legend at 0x7f9af6111130>


from risset import *

def do_risset_faster(filename, tune_length, freqs_per_note, sr):
    """
    Implement the faster version of Risset beats that aggregates
    duplicate frequencies into a sine and cosine term, and then
    uses the fast inverse Discrete Fourier Transform to reconstruct
    the time signal

    Parameters
    ----------
    filename: string
        Path to file with the tune
    tune_length: float
        Length, in seconds, of the tune
    freqs_per_note: int
        Number of frequencies to use for each note
    sr: int
        The sample rate of the entire piece
    """
    ps, times = load_tune(filename, tune_length)
    N = int(sr*tune_length)
    c = np.zeros(N) # Analagous to dictionaries
    s = np.zeros(N) # but using units of cycles per interval (instead of hz, cycles/sec)
    for p, shift in zip(ps, times):
        fhz = 440*(2**(p/12)) # In cycles/second (hz)
        k_center = int(fhz*N/sr)
        # Goal with risset beats to create a beat frequency that happens once
        # over our tune interval
        # Frequency k_center + 1 cycles/interval, difference is 1 cycle/interval
        for dk in range(-freqs_per_note//2, freqs_per_note//2+1):
            k = k_center + dk
            f = k*sr/N
            c[k] += np.cos(2*np.pi*f*shift)
            s[k] += np.sin(2*np.pi*f*shift)
    
    ## Setting up format to use the "Fast Fourier Transform," which is an O(N log N) algorithm
    ## whereas the "brute force" version we implemented is O(N^2)
    M = N//2
    if N%2 == 1:
        M = (N-1)//2
    else:
        M = (N-2)//2
    c[-M::] = c[1:M+1][::-1]
    s[-M::] = -s[1:M+1][::-1]
    return np.real(np.fft.ifft(c - 1j*s))


import time
tune_length = 70
freqs_per_note = 500
sr = 44100
tic = time.time()
x = do_risset_faster("Tunes/wanna.txt", tune_length, freqs_per_note, sr)
print("Elapsed Time: {:.3f}".format(time.time()-tic))
ipd.Audio(x, rate=sr)

Elapsed Time: 0.241


N = 44100*70
# N^2 / N log2(N),   N/log2(N)
N/np.log2(N)

143196.6024188806

The Inverse Real Discrete Fourier Transform¶

Forward Discrete Fourier Transform¶

$c[k] = \sum_{n=0}^{N-1} x[n] \cos(2 \pi k n / N) / (N/2)$¶

$s[k] = \sum_{n=0}^{N-1} x[n] \sin(2 \pi k n / N) / (N/2)$¶

Inverse Discrete Fourier Transform¶

$ x[n] = \sum_{k = 0}^{N/2} c[k] \cos(2 \pi k n / N) + s[k] \sin(2 \pi k n / N)$¶

Fast(er) Risset Beats¶